### Please also see README for required R packages
packages <- c(
  "tidyverse",   # includes dplyr, ggplot2, readr, stringr, forcats, purrr etc.
  "dplyr",       # for select
  "firatheme",   # for theme_fira()
  "ggplot2",     # for plotting
  "broom",       # for tidy model output (optional)
  "MASS",        # for some modeling functions (used indirectly)
  "scales",      # for axis formatting
  "stringr",     # string manipulation
  "forcats",     # factor handling
  "gridExtra",   # combining plots (optional)
  "cowplot",     # combining plots (optional)
  "haven",       # for read_sav() to load SPSS data
  "ggalluvial",  # for geom_alluvium() used in multiple figures
  "readr"        # for read_csv() used in Figure 2
)

lapply(packages, library, character.only = TRUE)


### Figure 1 - left panel
### load data
load("data/fig1_left.RData")

ggplot(processed_data, aes(x = Date)) +
  geom_rect(aes(xmin = as.Date("2022-03-06"), xmax = as.Date("2022-06-01"), ymin = 0, ymax = 60), fill = "lightgray", alpha = 0.5) +
  geom_vline(xintercept = as.numeric(as.Date("2022-04-22")), 
             color = "darkgray", linetype = "dashed", linewidth = 1) +
  geom_vline(xintercept = as.numeric(as.Date("2022-05-20")), 
             color = "darkgray", linetype = "dashed", linewidth = 1) +
  geom_bar(stat = "count", color = "gray30") +
  xlab("Date") +
  ylab("Number of News Articles") +
  ggtitle("") + theme_fira() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))



### Figure 1 - right panel
### load data
load("data/fig1_right.RData")

ggplot(processed_data, aes(x = Date)) +
  geom_rect(aes(xmin = as.Date("2022-03-06"), xmax = as.Date("2022-06-01"), ymin = 0, ymax = 60), fill = "lightgray", alpha = 0.5) +
  geom_vline(xintercept = as.numeric(as.Date("2022-04-22")), 
             color = "darkgray", linetype = "dashed", linewidth = 1) +
  geom_vline(xintercept = as.numeric(as.Date("2022-05-20")), 
             color = "darkgray", linetype = "dashed", linewidth = 1) +
  geom_bar(stat = "count", color = "gray30") +
  xlab("Date") +
  ylab("Number of TV Newscasts") +
  ggtitle("") + theme_fira() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

### Figure 2
### load data
topic_distr <- read_csv("data/topic_distr.csv")

df <- topic_distr %>%
  mutate(topic = pmap_chr(dplyr::select(., starts_with("t_")), 
                          ~which.max(c(...)) %>%
                            paste0("t_", .))) %>%
  mutate(topic_index = readr::parse_number(topic))

topic_distribution <- df %>%
  count(topic) %>%
  mutate(percentage = n / sum(n) * 100)

topic_distribution <- topic_distribution %>%
  mutate(topic = reorder(topic, -percentage))

ggplot(topic_distribution, aes(x = topic, y = percentage)) +
  geom_bar(stat = "identity", show.legend = FALSE, color = "gray20", fill = "gray80") +
  ylim(0,12) + theme_fira() +
  labs(x = "Topic", y = "Percentage (%)", title = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")

### Figure 3
### load data
load("data/panel_data.RData")

viden <- panel_pre_post %>%
  group_by(Q25_pre, Q25_post_post) %>%
  summarise(cnt = n()) %>%
  mutate(freq = round(cnt / sum(cnt), 3)) %>% 
  arrange(desc(freq))

viden$Q25_pre <- haven::as_factor(viden$Q25_pre)
viden$Q25_post_post <- haven::as_factor(viden$Q25_post_post)

viden$Q25_pre <- recode(viden$Q25_pre, 
                        "Helt enig" = "Completely agree", "Delvis enig" = "Partially agree", "Hverken enig eller uenig" = "Neither agree nor disagree",
                        "Delvis uenig" = "Partially disagree", "Helt uenig" = "Completely disagree")

viden$Q25_post_post <- recode(viden$Q25_post_post, 
                              "Helt enig" = "Completely agree", "Delvis enig" = "Partially agree", "Hverken enig eller uenig" = "Neither agree nor disagree",
                              "Delvis uenig" = "Partially disagree", "Helt uenig" = "Completely disagree")

ggplot(viden,
       aes(y = cnt,
           axis1 = Q25_pre, axis2 = Q25_post_post, label = Q25_post_post)) +
  geom_alluvium(aes(fill = Q25_post_post), curve_type = "sigmoid") +
  geom_stratum(aes(fill = after_stat(stratum)), size = 0.15) +
  geom_text(stat = "stratum", colour = "white", fontface = "bold", size = 1.25,
            aes(label = after_stat(stratum))) +
  scale_fill_manual(values = c("#045F5F",  "#E55451",  "#ECBA2A", "#3975D3", "pink")) +
  scale_x_discrete(limits = c("Pre-Election Survey", "Post-Election Survey")) + theme_fira() + 
  theme(legend.position = "none", axis.text.y=element_blank(), panel.background=element_blank(),
        panel.border=element_blank(),panel.grid.major=element_blank(), 
        panel.grid.minor=element_blank(),plot.background=element_blank(), 
        axis.line=element_blank(), plot.title = element_text(size=9, hjust = 0, face= "italic")) + ylab("") + 
  labs(title = "I have/had sufficient knowledge to assess what I should \nvote for in the referendum on the defense opt-out.")


### Figure 4
load("data/panel_data.RData")

df <- panel_pre_post %>%
  group_by(Q5_pre, Q5_post_post) %>%
  summarise(cnt = n()) %>%
  mutate(freq = round(cnt / sum(cnt), 3)) %>% 
  arrange(desc(freq))

df$Q5_pre <- haven::as_factor(df$Q5_pre)
df$Q5_post_post <- haven::as_factor(df$Q5_post_post)

df$Q5_post_post <- recode(df$Q5_post_post, 
                          'Ved ikke / vil ikke svare' = 'Ved ikke / vil ikke svare',
                          'Stemte ikke' = 'Stemte ikke / blankt',
                          'Stemte blankt' = 'Stemte ikke / blankt')

df$Q5_pre <- recode(df$Q5_pre, 
                    'Vil ikke stemme / vil stemme blankt' = 'Stemmer ikke / blankt')

df$Q5_post_post <- recode(df$Q5_post_post, 
                          "For (Ja)" = "Yes",                  "Imod (Nej)" = "No",               
                          "Stemte ikke / blankt" = "abstain",      "Ved ikke / vil ikke svare" = "don't know")

df$Q5_pre <- recode(df$Q5_pre, 
                    "For" = "Yes", "Imod" = "No", "Stemmer ikke / blankt" = "abstain",
                    "Ved ikke" = "don't know"     )

ggplot(df,
       aes(y = cnt,
           axis1 = Q5_pre, axis2 = Q5_post_post, label = Q5_post_post)) +
  geom_alluvium(aes(fill = Q5_post_post), curve_type = "sigmoid") +
  geom_stratum(aes(fill = after_stat(stratum)), size = 0.15) +
  geom_text(stat = "stratum", colour = "white", fontface = "bold", size = 1.25,
            aes(label = after_stat(stratum))) +
  scale_fill_manual(values = c("#045F5F", "#E55451",  "#ECBA2A", "gray80")) +
  scale_x_discrete(limits = c("Pre-Election Survey", "Post-Election Survey")) + theme_fira() + 
  theme(legend.position = "none", axis.text.y=element_blank(), panel.background=element_blank(),
        panel.border=element_blank(),panel.grid.major=element_blank(), 
        panel.grid.minor=element_blank(),plot.background=element_blank(), 
        axis.line=element_blank(), plot.title = element_text(size=9, hjust = 0, face= "italic")) + ylab("") + 
  labs(title = "Do you vote/Did you vote for (yes) or against (no) the abolition of the Danish \nopt-out from the EU's common defense and security policy \nat the referendum on 1 June?")

### Figure 5 - upper panel 
load("data/panel_data.RData")

eu_position <- panel_pre_post |>
  mutate(
    Q15_pre_recoded = case_when(
      Q15_pre %in% c("1", "2") ~ "More \nadvantages",
      Q15_pre %in% c("3", "4") ~ "More \ndisadvantages",
      TRUE ~ NA_character_
    ),
    Q15_post_recoded = case_when(
      Q15_post %in% c("1", "2") ~ "More \nadvantages",
      Q15_post %in% c("3", "4") ~ "More \ndisadvantages",
      TRUE ~ NA_character_
    )
  ) |>
  filter(!is.na(Q15_pre_recoded), !is.na(Q15_post_recoded)) |>
  count(Q15_pre_recoded, Q15_post_recoded, name = "cnt") |>
  mutate(freq = round(cnt / sum(cnt), 3)) |>
  arrange(desc(freq))

ggplot(eu_position,
       aes(y = cnt,
           axis1 = Q15_pre_recoded, axis2 = Q15_post_recoded, label = Q15_post_recoded)) +
  geom_alluvium(aes(fill = Q15_post_recoded), curve_type = "sigmoid") +
  geom_stratum(aes(fill = after_stat(stratum)), size = 0.15) +
  geom_text(stat = "stratum", colour = "white", fontface = "bold", size = 2,
            aes(label = after_stat(stratum))) +
  scale_fill_manual(values = c("#045F5F", "#E55451")) +
  scale_x_discrete(limits = c("Pre-Election Survey", "Post-Election Survey")) + theme_fira() + 
  theme(legend.position = "none", axis.text.y=element_blank(), panel.background=element_blank(),
        panel.border=element_blank(),panel.grid.major=element_blank(), 
        panel.grid.minor=element_blank(),plot.background=element_blank(), 
        axis.line=element_blank(), plot.title = element_text(size=9, hjust = 0, face= "italic")) + ylab("") + 
  labs(title = "All things considered, do you think Denmark has had more advantages \nor more disadvantages from its EU membership?")


### Figure 5 - lower panel 
load("data/panel_data.RData")

eu_vote_assoc <- panel_pre_post |>
  mutate(
    # Recode EU attitudes
    EU_pre = case_when(
      Q15_pre %in% c("1", "2") ~ "More \nadvantages",
      Q15_pre %in% c("3", "4") ~ "More \ndisadvantages",
      TRUE ~ NA_character_
    ),
    EU_post = case_when(
      Q15_post %in% c("1", "2") ~ "More \nadvantages",
      Q15_post %in% c("3", "4") ~ "More \ndisadvantages",
      TRUE ~ NA_character_
    ),
    # Recode all voting behavior categories with clear labels
    Vote_pre = recode(as.character(Q5_pre),
                      "1" = "Yes (For)", 
                      "2" = "No (Against)",
                      "3" = "Abstained/Blank",
                      "4" = "Don't know"),
    Vote_post = recode(as.character(Q5_post_post),
                       "1" = "Yes (For)", 
                       "2" = "No (Against)",
                       "3" = "Abstained/Blank",
                       "4" = "Don't know",
                       "5" = NA_character_)  # Drop category 5
  ) |>
  # Remove rows where Vote_post was category 5
  filter(!is.na(EU_pre), !is.na(EU_post), !is.na(Vote_pre), !is.na(Vote_post)) |>
  # Gather data for plotting
  pivot_longer(cols = c(EU_pre, EU_post), names_to = "Survey_Time", values_to = "EU_Attitude") |>
  pivot_longer(cols = c(Vote_pre, Vote_post), names_to = "Vote_Time", values_to = "Vote") |>
  filter((Survey_Time == "EU_pre" & Vote_Time == "Vote_pre") | 
           (Survey_Time == "EU_post" & Vote_Time == "Vote_post")) |>
  count(Survey_Time, EU_Attitude, Vote, name = "cnt") |>
  group_by(Survey_Time, EU_Attitude) |>
  mutate(percent = round(100 * cnt / sum(cnt), 1)) |>
  # Ensure correct facet order
  mutate(Survey_Time = factor(Survey_Time, levels = c("EU_pre", "EU_post")))


ggplot(eu_vote_assoc, aes(x = EU_Attitude, y = percent, fill = Vote)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  
  # Add percentage labels on the bars
  geom_text(aes(label = paste0(percent, "%")), 
            position = position_dodge(width = 0.7), 
            vjust = -0.5, size = 3) +
  
  # Facet to show Pre- and Post-Referendum responses
  facet_wrap(~ Survey_Time, 
             labeller = as_labeller(c("EU_pre" = "Pre-Referendum", "EU_post" = "Post-Referendum"))) +
  
  # Define colors for all voting categories
  scale_fill_manual(
    values = c(
      "Yes (For)" = "#045F5F",      # Dark Green
      "No (Against)" = "#E55451",   # Red
      "Abstained/Blank" = "#ECBA2A",# Yellow
      "Don't know" = "gray80"       # Gray
    )
  ) +
  
  theme_fira() +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    axis.title.x = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.background = element_blank(),
    plot.title = element_text(size = 16, hjust = 0, face = "italic")
  ) +
  
  ylab("Percentage of Respondents") +
  labs(title = "Association Between EU Attitudes and Voting Behavior Before and After the Referendum")

### Figure 6 - upper left panel
### load data
df <- read_sav("data/Survey_ex_1_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0

# factor and re-level treatment variable
df$treatment <- as.factor(df$QGroup_Rand)
table(df$treatment)
df$treatment <- relevel(df$treatment, ref = 5)
table(df$treatment)

#treatment collapsed
df <- df %>%
  mutate(
    treatment_collapsed = case_when(
      treatment %in% c(1, 2) ~ "Yes",
      treatment %in% c(3, 4) ~ "No",
      treatment == 5 ~ "Control",
    )) 


# Base model with collapsed treatment
model_collapsed <- lm(vote_yes ~ treatment_collapsed, data = df, weights = weight1)
summary(model_collapsed)
# Extract confidence intervals
confint95_collapsed <- confint(model_collapsed, level = 0.95)
confint90_collapsed <- confint(model_collapsed, level = 0.90)

# Extract estimates and confidence intervals (excluding intercept)
estimates_collapsed <- model_collapsed$coefficients[-1] * 100
lower95_collapsed <- confint95_collapsed[-1, "2.5 %"] * 100
upper95_collapsed <- confint95_collapsed[-1, "97.5 %"] * 100
lower90_collapsed <- confint90_collapsed[-1, "5 %"] * 100
upper90_collapsed <- confint90_collapsed[-1, "95 %"] * 100

# Define labels
labels_collapsed <- c("No \n(vs. Control)", "Yes \n(vs. Control)")

# Ensure correct data frame construction
plotdf_collapsed <- data.frame(
  estimates = estimates_collapsed, 
  lower95 = lower95_collapsed, 
  upper95 = upper95_collapsed, 
  lower90 = lower90_collapsed, 
  upper90 = upper90_collapsed, 
  labels = labels_collapsed
)

# Main plot
ggplot(plotdf_collapsed, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-12, 12) +
  scale_color_manual(values=c("gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "") +
  ggtitle("Start of campaign \n(Collapsed)") +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()


### Figure 6 - upper right panel 
### load data
df <- read_sav("data/Survey_ex_2_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0


# factor and re-level treatment variable
df$treatment <- as.factor(df$QGroup_Rand)
table(df$treatment)
df$treatment <- relevel(df$treatment, ref = 5)
table(df$treatment)


#treatment collapsed
df <- df %>%
  mutate(
    treatment_collapsed = case_when(
      treatment %in% c(1, 2) ~ "Yes",
      treatment %in% c(3, 4) ~ "No",
      treatment == 5 ~ "Control",
    )) 

table(df$treatment_collapsed)

# Base model with collapsed treatment
model_collapsed <- lm(vote_yes ~ treatment_collapsed, data = df, weights = weight1)
summary(model_collapsed)
# Extract confidence intervals
confint95_collapsed <- confint(model_collapsed, level = 0.95)
confint90_collapsed <- confint(model_collapsed, level = 0.90)

# Extract estimates and confidence intervals (excluding intercept)
estimates_collapsed <- model_collapsed$coefficients[-1] * 100
lower95_collapsed <- confint95_collapsed[-1, "2.5 %"] * 100
upper95_collapsed <- confint95_collapsed[-1, "97.5 %"] * 100
lower90_collapsed <- confint90_collapsed[-1, "5 %"] * 100
upper90_collapsed <- confint90_collapsed[-1, "95 %"] * 100

# Define labels
labels_collapsed <- c("No \n(vs. Control)", "Yes \n(vs. Control)")

# Ensure correct data frame construction
plotdf_collapsed <- data.frame(
  estimates = estimates_collapsed, 
  lower95 = lower95_collapsed, 
  upper95 = upper95_collapsed, 
  lower90 = lower90_collapsed, 
  upper90 = upper90_collapsed, 
  labels = labels_collapsed
)

# Main plot
ggplot(plotdf_collapsed, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-12, 12) +
  scale_color_manual(values=c("gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "") +
  ggtitle("End of campaign \n(Collapsed)") +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()


### Figure 6 - lower left panel  
### load data
df <- read_sav("data/Survey_ex_1_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0

# factor and re-level treatment variable
df$treatment <- as.factor(df$QGroup_Rand)
table(df$treatment)
df$treatment <- relevel(df$treatment, ref = 5)
table(df$treatment)

# base model
summary(lm(vote_yes ~ treatment, data = df, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = df, weights = weight1)
confint95 <- as.data.frame(confint(model1, level = 0.95))
confint90 <- as.data.frame(confint(model1, level = 0.90))

  
df %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# based on base model (see above)
estimates <- model1$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf1 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf1, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-12, 12) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "") +
  ggtitle("Start of campaign") +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()

### Figure 6 - lower right panel 
### load data
df <- read_sav("data/Survey_ex_2_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0


# factor and re-level treatment variable
df$treatment <- as.factor(df$QGroup_Rand)
table(df$treatment)
df$treatment <- relevel(df$treatment, ref = 5)
table(df$treatment)

# base model
summary(lm(vote_yes ~ treatment, data = df, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = df, weights = weight1)
confint95 <- as.data.frame(confint(model1, level = 0.95))
confint90 <- as.data.frame(confint(model1, level = 0.90))


df %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# based on base model (see above)
estimates <- model1$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf1 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf1, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-12, 12) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "") +
  ggtitle("End of campaign") + coord_flip() +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank())

### Figure 7 - left panel 
### load data
df <- read_sav("data/Survey_ex_1_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0

# factor and re-level treatment variable
df$treatment <- as.factor(df$QGroup_Rand)
table(df$treatment)
df$treatment <- relevel(df$treatment, ref = 5)
table(df$treatment)

#subset data by excluding don't know answers
df_subset <- subset(df, Q8 == 1 | Q8 == 2 | Q8 == 3)
df_subset$vote_yes <- NA
df_subset$vote_yes[df_subset$Q8 == 1] <- 1
df_subset$vote_yes[df_subset$Q8 != 1] <- 0
table(df_subset$vote_yes)
table(df_subset$Q8)

# factor and re-level treatment variable
df_subset$treatment <- as.factor(df_subset$QGroup_Rand)
table(df_subset$treatment)
df_subset$treatment <- relevel(df_subset$treatment, ref = 5)

# base model
summary(lm(vote_yes ~ treatment, data = df_subset, weights = weight1))
model3 <- lm(vote_yes ~ treatment, data = df_subset, weights = weight1)
confint95 <- as.data.frame(confint(model3, level = 0.95))
confint90 <- as.data.frame(confint(model3, level = 0.90))

# based on base model (see above)
estimates <- model3$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf3 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf3, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-16, 16) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "") +
  ggtitle("Start of campaign") + coord_flip() +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) 


### Figure 7 - right panel 
### load data
df <- read_sav("data/Survey_ex_2_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0

# factor and re-level treatment variable
df$treatment <- as.factor(df$QGroup_Rand)
df$treatment <- relevel(df$treatment, ref = 5)

#subset data by excluding don't know answers
df_subset <- subset(df, Q8 == 1 | Q8 == 2 | Q8 == 3)
df_subset$vote_yes <- NA
df_subset$vote_yes[df_subset$Q8 == 1] <- 1
df_subset$vote_yes[df_subset$Q8 != 1] <- 0
table(df_subset$vote_yes)
table(df_subset$Q8)

# factor and re-level treatment variable
df_subset$treatment <- as.factor(df_subset$QGroup_Rand)
table(df_subset$treatment)
df_subset$treatment <- relevel(df_subset$treatment, ref = 5)

df_subset %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# base model
summary(lm(vote_yes ~ treatment, data = df_subset, weights = weight1))
model3 <- lm(vote_yes ~ treatment, data = df_subset, weights = weight1)
confint95 <- as.data.frame(confint(model3, level = 0.95))
confint90 <- as.data.frame(confint(model3, level = 0.90))

# based on base model (see above)
estimates <- model3$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf3 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf3, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-16, 16) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "") +
  ggtitle("End of campaign") + coord_flip() +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank())  

### Table 2
### load data
df1 <- read_sav("data/Survey_ex_1_weigthed.sav")

# make binary outcome var
df1$vote_yes <- NA
df1$vote_yes[df1$Q8 == 1] <- 1
df1$vote_yes[df1$Q8 != 1] <- 0

# factor and re-level treatment variable
df1$treatment <- as.factor(df1$QGroup_Rand)
df1$treatment <- relevel(df1$treatment, ref = 5)

# base model
summary(lm(vote_yes ~ treatment, data = df1, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = df1, weights = weight1)

### load data
df2 <- read_sav("data/Survey_ex_2_weigthed.sav")

# make binary outcome var
df2$vote_yes <- NA
df2$vote_yes[df2$Q8 == 1] <- 1
df2$vote_yes[df2$Q8 != 1] <- 0

# factor and re-level treatment variable
df2$treatment <- as.factor(df2$QGroup_Rand)
df2$treatment <- relevel(df2$treatment, ref = 5)

# base model
summary(lm(vote_yes ~ treatment, data = df2, weights = weight1))
model2 <- lm(vote_yes ~ treatment, data = df2, weights = weight1)

extract_z_p <- function(model1, model2) {
  # Initialize vectors
  positions <- 2:5
  results <- data.frame(
    Position = positions,
    Coef_Model1 = NA,
    Coef_Model2 = NA,
    Z = NA,
    P = NA
  )
  
  for (i in positions) {
    coef1 <- coef(model1)[i]
    coef2 <- coef(model2)[i]
    se1 <- summary(model1)$coefficients[i, 2]
    se2 <- summary(model2)$coefficients[i, 2]
    
    z <- (coef1 - coef2) / sqrt(se1^2 + se2^2)
    p <- 2 * (1 - pnorm(abs(z)))
    
    results[results$Position == i, "Coef_Model1"] <- coef1
    results[results$Position == i, "Coef_Model2"] <- coef2
    results[results$Position == i, "Z"] <- z
    results[results$Position == i, "P"] <- p
  }
  
  return(results)
}

extract_z_p(model1, model2)

### Figure 8 - left panel 
### load data
df <- read_sav("data/Survey_ex_1_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df1$Q8 == 1] <- 1
df$vote_yes[df1$Q8 != 1] <- 0

# factor and re-level treatment variable
df$treatment <- as.factor(df1$QGroup_Rand)
df$treatment <- relevel(df1$treatment, ref = 5)

# subsetting for lower education
loweducation <- df %>% 
  filter(Q4 < 6)

# factor and re-level treatment variable
loweducation$treatment <- as.factor(loweducation$QGroup_Rand)
table(loweducation$treatment)
loweducation$treatment <- relevel(loweducation$treatment, ref = 5)
table(loweducation$treatment)

# base model
summary(lm(vote_yes ~ treatment, data = loweducation, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = loweducation, weights = weight1)
confint95 <- as.data.frame(confint(model1, level = 0.95))
confint90 <- as.data.frame(confint(model1, level = 0.90))

loweducation %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# based on base model (see above)
estimates <- model1$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf1 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf1, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-20, 20) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="Less educated", caption = "") +
  theme(axis.title = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray20", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()

### Figure 8 - right panel 
### load data
df <- read_sav("data/Survey_ex_1_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0


# subsetting for more educated
education <- df %>% 
  filter(Q4 > 5)

# factor and re-level treatment variable
education$treatment <- as.factor(education$QGroup_Rand)
education$treatment <- relevel(education$treatment, ref = 5)

# base model
summary(lm(vote_yes ~ treatment, data = education, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = education, weights = weight1)
confint95 <- as.data.frame(confint(model1, level = 0.95))
confint90 <- as.data.frame(confint(model1, level = 0.90))

education %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# based on base model (see above)
estimates <- model1$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf1 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf1, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-20, 20) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="More educated", caption = "") +
  theme(axis.title = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()

### Figure 9 - left panel 
### load data
df <- read_sav("data/Survey_ex_2_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0

loweducation <- df %>% 
  filter(Q4 < 6)

# factor and re-level treatment variable
loweducation$treatment <- as.factor(loweducation$QGroup_Rand)
loweducation$treatment <- relevel(loweducation$treatment, ref = 5)

# base model
summary(lm(vote_yes ~ treatment, data = loweducation, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = loweducation, weights = weight1)
confint95 <- as.data.frame(confint(model1, level = 0.95))
confint90 <- as.data.frame(confint(model1, level = 0.90))

loweducation %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# based on base model (see above)
estimates <- model1$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf1 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf1, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-20, 20) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="Less educated", caption = "") +
  theme(axis.title = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray20", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()

### Figure 9 - left panel 
### load data
df <- read_sav("data/Survey_ex_2_weigthed.sav")

# make binary outcome var
df$vote_yes <- NA
df$vote_yes[df$Q8 == 1] <- 1
df$vote_yes[df$Q8 != 1] <- 0

# subsetting for more educated
education <- df %>% 
  filter(Q4 > 5)

# factor and re-level treatment variable
education$treatment <- as.factor(education$QGroup_Rand)
education$treatment <- relevel(education$treatment, ref = 5)

# base model
summary(lm(vote_yes ~ treatment, data = education, weights = weight1))
model1 <- lm(vote_yes ~ treatment, data = education, weights = weight1)
confint95 <- as.data.frame(confint(model1, level = 0.95))
confint90 <- as.data.frame(confint(model1, level = 0.90))

education %>%
  group_by(treatment) %>%
  count(vote_yes) %>%
  mutate(prop = prop.table(n))

# based on base model (see above)
estimates <- model1$coefficients[2:5]*100
lower95 <- confint95$`2.5 %`[2:5]*100
upper95 <- confint95$`97.5 %`[2:5]*100
lower90 <- confint90$`5 %`[2:5]*100
upper90 <- confint90$`95 %`[2:5]*100
labels <- c("Yes Gain \n(vs. control)", "Yes Loss \n(vs. control)", "No Loss \n(vs. control)", "No Gain \n(vs. control)")
plotdf1 <- data.frame(estimates, lower95, upper95, lower90, upper90, labels)

#main plot
ggplot(plotdf1, aes(x=labels, y=estimates, group=labels, color = labels)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                            color = "gray50", size=0.5) + ylim(-25, 25) +
  scale_color_manual(values=c("gray20", "gray20", "gray20", "gray20")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="More educated", caption = "") +
  theme(axis.title = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = c(0)), strip.text.x = element_blank()) + coord_flip()

